To see how we got here, check the notebooks with the loading the GFS atmospheric forecast data, the loading of the buoy data, and the evaluation of the slab-ocean model.
Details: https://github.com/willirath/nia-prediction-low-latitudes/
# parameters
# regional coverage
lat_min, lat_max = -35, 35
lon_min, lon_max = 0, 360
# good forecasts reach
good_forecast_days = 7 # 7 days
# additional buoy positions to be plotted
# The ones found in the downloaded data will be shown anyway.
added_buoy_positions = [
{"lat": 20.0, "lon": -38.0},
{"lat": 15.0, "lon": -38.0},
{"lat": 21.0, "lon": -23.0},
{"lat": 12.0, "lon": -23.0},
{"lat": -6.0, "lon": -10.0},
{"lat": -10.0 ,"lon": -10.0},
]
# data files
GFS_zarr_store = "tmp_GFS.zarr"
slab_zarr_store = "tmp_slab.zarr"
buoy_file_name = "tmp_buoy_data"
buoy_positions_file = "tmp_buoy_positions.csv"
mimoc_mld_file = "tmp_mimoc_mld.nc"
# dask specifics
dask_kwargs = {"n_workers": 1, "threads_per_worker": 2, "memory_limit": 6e9}
Before doing any calculations, we'll need to import a few modules. We'll also start a Dask cluster for parallel execution.
# dask
from dask.distributed import Client
# plotting
from bokeh.models.formatters import DatetimeTickFormatter
import cartopy.crs as ccrs
import cmocean
import geoviews as gv
import holoviews as hv
import hvplot.xarray, hvplot.pandas
# numerics
import numpy as np
import pandas as pd
import xarray as xr
# aux
from functools import reduce
from operator import add
# create Dask cluster
client = Client(**dask_kwargs)
client
buoy_positions = pd.read_csv(buoy_positions_file)
added_buoy_positions =pd.DataFrame.from_records(added_buoy_positions)
added_buoy_positions["lon"] = (360.0 + added_buoy_positions["lon"]) % 360.0
buoy_positions = buoy_positions.merge(added_buoy_positions, how="outer")
buoy_positions.head(3)
ds_GFS = xr.open_zarr(GFS_zarr_store)
ds_slab = xr.open_zarr(slab_zarr_store)
ds_mld = xr.open_dataset(mimoc_mld_file)
We'll need the time stamp of the start of the forecasting data.
start_of_forecast = (~ds_GFS["is_forecast"].astype(bool)).sum().compute().data
start_of_forecast = ds_GFS["time"].data[max(0, start_of_forecast-1)]
print(start_of_forecast)
good_forecast_time = np.timedelta64(good_forecast_days, "D")
ds_GFS = ds_GFS.sel(
lat=slice(lat_max, lat_min),
lon=slice(lon_min, lon_max),
)
ds_slab = ds_slab.sel(
lat=slice(lat_max, lat_min),
lon=slice(lon_min, lon_max),
)
ds_slab
buoy_positions = buoy_positions.where(
buoy_positions["lat"].apply(pd.Interval(lat_min, lat_max).__contains__)
& buoy_positions["lon"].apply(pd.Interval(lon_min, lon_max).__contains__)
).dropna()
years_data = list(ds_slab.time.groupby("time.year").groups.keys())
years_min = min(years_data)
years_max = min(years_data)
duration_years_data = years_max - years_min + 1
# pad by one
years_min -= 1
duration_years_data += 2
mld_time_coord_lower = xr.DataArray(
[
np.datetime64(f"{years_min + m // 12:04d}-{(m % 12) + 1:02d}-01")
for m in range(0, duration_years_data * 12)
],
dims=("time", )
)
mld_time_coord_upper = xr.DataArray(
[
np.datetime64(f"{years_min + m // 12:04d}-{(m % 12) + 1:02d}-01")
for m in range(1, duration_years_data * 12 + 1)
],
dims=("time", )
)
mld_time_coord = (
mld_time_coord_lower
+ (mld_time_coord_upper - mld_time_coord_lower) / 2.0
)
display(mld_time_coord)
mld_expand = xr.concat(
[ds_mld for n in range(duration_years_data)],
dim="month"
).rename({"month": "time"})
mld_expand.coords["time"] = mld_time_coord
mld_expand
mld_slab = mld_expand.interp_like(
ds_slab.coords["time"]
).sel(
lat=ds_slab.coords["lat"], lon=ds_slab.coords["lon"], method="nearest"
)
mld_slab = mld_slab.compute()
display(mld_slab)
ds_slab["u_slab"] *= ds_slab.attrs["slab_model_H"] / mld_slab.mixed_layer_depth
ds_slab["v_slab"] *= ds_slab.attrs["slab_model_H"] / mld_slab.mixed_layer_depth
ds_slab["umag_slab"] *= ds_slab.attrs["slab_model_H"] / mld_slab.mixed_layer_depth
We'll plot the time-maximum of near-inertial speed for the good forecast period and for the whole forecast period.
First, we construct the plot.
slab_umag_good_forecast_max = ds_slab["umag_slab"].sel(
time=slice(start_of_forecast, start_of_forecast + good_forecast_time)
).max("time")
slab_umag_whole_forecast_max = ds_slab["umag_slab"].sel(
time=slice(start_of_forecast, None)
).max("time")
near_inertial_max_plots = (
(
slab_umag_good_forecast_max.hvplot(
x="lon", y="lat", z="umag_slab",
clim=(0, 1.0),
cmap=cmocean.cm.speed,
frame_width=800,
hover=False,
geo=True, coastline=True,
crs=ccrs.PlateCarree(), projection=ccrs.PlateCarree(),
title="Near-inertial speed max [m/s] and mixed-layer-depth [m], good forecast period"
)
+ slab_umag_whole_forecast_max.hvplot(
x="lon", y="lat", z="umag_slab",
clim=(0, 1.0),
cmap=cmocean.cm.speed,
frame_width=800,
hover=False,
geo=True, coastline=True,
crs=ccrs.PlateCarree(), projection=ccrs.PlateCarree(),
title="Near-inertial speed max [m/s] and mixed-layer-depth [m], whole forecast period"
)
) * mld_slab.mean("time").mixed_layer_depth.hvplot.contour(
x="lon", y="lat", geo=True, cmap="gray", hover=True,
levels=list(range(0, 50, 10)) + list(range(60, 120, 20)), line_width=1.5, alpha=0.5
) * buoy_positions.hvplot.points(
y="lat", x="lon", geo=True, coastline=True,
marker='diamond',
fill_color=None, line_color="black",
line_width=2, size=70,
) * gv.feature.grid()
).cols(1)
display(near_inertial_max_plots)
time_series_plots = []
time_formatter = DatetimeTickFormatter(
months='%b %Y', days='%b %d'
)
forecast_spans = (
hv.VSpan(
start_of_forecast, start_of_forecast + good_forecast_time
).opts(padding=0, color='lightgray')
* hv.VSpan(
start_of_forecast + good_forecast_time, None
).opts(padding=0, color='pink')
)
for lat, lon in zip(buoy_positions["lat"], buoy_positions["lon"]):
local_mld = mld_slab.mixed_layer_depth.sel(lat=lat, lon=lon, method='nearest').mean('time').data
name = f"{lat}N {lon}E, MLD={local_mld:.0f}m"
buoy_ds = ds_slab.sel(lat=lat, lon=lon, method="nearest")
buoy_ds["U20"] = ds_GFS["U20"].sel(lat=lat, lon=lon, method="nearest")
buoy_ds["V20"] = ds_GFS["V20"].sel(lat=lat, lon=lon, method="nearest")
if (buoy_ds["umag_slab"].max("time").isnull().data.compute()):
continue
time_series_plots.append(
(
(
forecast_spans.redim.label(y="u_slab")
* buoy_ds["u_slab"].hvplot.line(label="zonal near-inertial current")
* buoy_ds["v_slab"].hvplot.line(label="meridional near-inertial current")
* buoy_ds["umag_slab"].hvplot.line(label="near-inertial speed")
).options(
width=800, height=160, show_grid=True,
xaxis=None,
legend_cols=False, legend_position='right',
ylabel="current [m/s]", title=name
)
+ (
forecast_spans.redim.label(y="U20")
* buoy_ds["U20"].hvplot.line(label="zonal wind (20m)")
* buoy_ds["V20"].hvplot.line(label="meridional wind (20m)")
).options(
width=800, height=160, show_grid=True,
xformatter=time_formatter,
legend_cols=False, legend_position='right',
ylabel="wind [m/s]", xlabel=""
)
)
)
time_series_plots = reduce(add, time_series_plots).cols(1)
display(time_series_plots)
To get a feeling for the atmospheric conditions, we'll plot sea-level pressure anomalies every 12 hours for 3 days before and throughout the whole forecast period.
Anomalies are calculated relative to the whole data period (usually 30+14 days).
SLP = ds_GFS["SLP"].compute()
SLP_mean = SLP.mean("time")
SLP_anomaly = (SLP - SLP_mean)
plot_every = np.timedelta64(12, "h")
max_iter = ((SLP_anomaly.coords["time"].max("time") - start_of_forecast) / plot_every).item() // 1 + 1
plot_times = [
(start_of_forecast + n * plot_every)
for n in range(-6, int(max_iter))
]
plots = []
for plot_time in plot_times:
title = f"SLP anomaly [hPa], {pd.Timestamp(plot_time).strftime('%Y-%m-%d %H:%M:%S UTC')}"
if plot_time > start_of_forecast:
title += f"\t(forecast + {(plot_time - start_of_forecast) / np.timedelta64(1, 'h')}h)"
try:
plots.append(
(
SLP_anomaly.sel(time=plot_time, method="nearest").compute().hvplot(
clim=(-10, 10),
cmap=cmocean.cm.delta,
frame_width=800,
geo=True, coastline=True,
crs=ccrs.PlateCarree(), projection=ccrs.PlateCarree(),
hover=False
)
* buoy_positions.hvplot.points(
y="lat", x="lon", geo=True, coastline=True,
marker='circle',
fill_color=None, line_color="black",
line_width=1.5, size=60,
)
* gv.feature.grid()
).opts(
title=title,
show_grid=True
)
)
except Exception as e:
print(f"for {plot_time} I got: {e}")
slp_plot = reduce(add, plots).cols(1)
display(slp_plot)